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Abstract 

We  present  spectral  element  (SE)  and  discontinuous  Galerkin  (DG)  solutions  of  the  Euler  and  compressible  Navier- 
Stokes  (NS)  equations  for  stratified  fluid  flow  which  are  of  importance  in  nonhydrostatic  mesoscale  atmospheric  modeling. 
We  study  three  different  forms  of  the  governing  equations  using  seven  test  cases.  Three  test  cases  involve  flow  over  moun¬ 
tains  which  require  the  implementation  of  non-reflecting  boundary  conditions,  while  one  test  requires  viscous  terms  (den¬ 
sity  current).  Including  viscous  stresses  into  finite  difference,  finite  element,  or  spectral  element  models  poses  no  additional 
challenges;  however,  including  these  terms  to  either  finite  volume  or  discontinuous  Galerkin  models  requires  the  introduc¬ 
tion  of  additional  machinery  because  these  methods  were  originally  designed  for  first-order  operators.  We  use  the  local 
discontinuous  Galerkin  method  to  overcome  this  obstacle.  The  seven  test  cases  show  that  all  of  our  models  yield  good 
results.  The  main  conclusion  is  that  equation  set  1  (non-conservation  form)  does  not  perform  as  well  as  sets  2  and  3  (con¬ 
servation  forms).  For  the  density  current  (viscous),  the  SE  and  DG  models  using  set  3  (mass  and  total  energy)  give  less 
dissipative  results  than  the  other  equation  sets;  based  on  these  results  we  recommend  set  3  for  the  development  of  future 
multiscale  research  codes.  In  addition,  the  fact  that  set  3  conserves  both  mass  and  energy  up  to  machine  precision  motives 
us  to  pursue  this  equation  set  for  the  development  of  future  mesoscale  models.  For  the  bubble  and  mountain  tests,  the  DG 
models  performed  better.  Based  on  these  results  and  due  to  its  conservation  properties  we  recommend  the  DG  method.  In 
the  worst  case  scenario,  the  DG  models  are  50%  slower  than  the  non-conservative  SE  models.  In  the  best  case  scenario,  the 
DG  models  are  just  as  efficient  as  the  conservative  SE  models. 
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1.  Introduction 

The  recent  paradigm  shift  in  high-performance  computing  (HPC)  is  forcing  many  numerical  weather  pre¬ 
diction  (NWP)  operational  centers  to  rethink  the  numerical  methods  that  their  models  are  based  on.  For 
example,  the  current  trend  in  distributed-memory  computing  has  moved  toward  clusters  based  on  hundreds 
of  thousands  of  cheap,  commodity-based  processors;  the  top  three  fastest  computers  in  the  world  in  2007  have 
212,000  (Lawrence  Livermore  National  Laboratory,  USA),  65,000  (Forschungszentrum  Juelich,  Germany), 
and  14,000  (New  Mexico  Computing  Applications  Center,  USA).  It  is  expected  that  clusters  comprised  of  mil¬ 
lions  of  processors  will  appear  very  soon.  In  order  to  take  full  advantage  of  computers  with  such  high  proces¬ 
sor  counts  requires  exploring  numerical  methods  that  are  local  in  nature,  have  a  large  on-processor  operation 
count,  and  a  small  communication  footprint.  Local  high-order  methods  like  the  spectral  element  and  discon¬ 
tinuous  Galerkin  methods  have  all  of  these  properties  and  for  this  reason  they  have  been  successfully  applied 
to  a  variety  of  problems. 

Spectral  element  (SE)  methods  combine  the  local  domain  decomposition  property  of  finite  element  (FE) 
methods  with  the  high-order  accuracy  and  weak  numerical  dispersion  of  spectral  methods.  SE  methods  have 
shown  promise  in  many  areas  of  the  geosciences  including:  seismic  wave  propagation  [33],  deep  Earth  flows 
[13],  climate  [53,14],  ocean  [28,38],  and  numerical  weather  prediction  [21,22].  These  methods  are  high-order 
FE  methods  where  the  interpolation  and  integration  points  are  chosen  to  be  the  Legendre-Gauss-Lobatto 
points. 

In  contrast,  discontinuous  Galerkin  (DG)  methods  combine  the  local  domain  decomposition  property  of 
FE  methods,  the  high-order  accuracy  and  weak  numerical  dispersion  of  spectral  methods,  and  the  conserva¬ 
tion  properties  of  finite  volume  (FV)  methods;  in  essence,  DG  methods  are  the  high-order  generalization  of 
FV  methods.  There  are  two  distinct  types  of  DG  methods:  nodal  (see  [20,23])  and  modal  (see  [12,59]),  but 
in  the  current  study  we  only  consider  the  nodal  approach  introduced  in  [20]  which  uses  the  same  machinery 
developed  for  SE  methods  such  as  quadrilateral  grids,  tensor  product  basis  functions,  and  Legendre-Gauss- 
Lobatto  grid  points.  It  has  been  only  very  recent  (since  2000)  that  the  DG  method  first  appeared  in  geophys¬ 
ical  fluid  dynamics  (GFD)  applications.  However,  implementations  of  the  DG  method  in  GFD  have  remained 
primarily  restricted  to  shallow  water  flow  (see  [44,35,2,20,10,12,37,40,34,23,24]).  To  date,  there  has  been  no 
published  work  on  either  SE  or  DG  methods  for  nonhydrostatic  mesoscale  atmospheric  applications. 

However,  doing  something  for  the  first  time  is  not  a  sufficient  reason  for  developing  a  new  model  -  the  new 
model  must  have  attractive  properties  not  offered  by  existing  models.  The  high-order  accuracy,  geometric  flex¬ 
ibility  to  use  any  grid,  and  the  scalability  of  SE  and  DG  methods  on  large  processor  count  computers  are  suf¬ 
ficient  reasons  for  exploring  this  new  class  of  models. 

Almost  all  nonhydrostatic  mesoscale  models  currently  in  existence  are  based  on  the  finite  difference  (FD) 
method;  examples  include  the  following  list  of  models:  [4,5,11,17,25,26,29-32,39,42,43,46,48,57,58],  and  [60]. 
Included  in  this  list  are  models  such  as  ARPS  (University  of  Oklahoma),  COAMPS  (US  Navy),  LM  (German 
Weather  Service),  MC2  (Environment  Canada),  MM5  (Penn  State/NCAR),  NMM  (National  Center  for 
Environmental  Prediction),  and  WRF  (NCAR).  The  only  models  in  the  literature  not  based  on  the  FD 
method  are  the  FV  models  found  in  [7]  and  [1],  and  the  SE  and  DG  models  presented  in  our  paper.  One 
of  the  biggest  advantages  that  SE  and  DG  methods  have  over  the  FD  method  is  that  no  terrain  following 
coordinates  of  the  type  presented  in  [16]  need  to  be  included  in  the  governing  equations.  Of  course,  the  orog¬ 
raphy  has  to  be  accounted  for  in  some  manner  but  element-based  Galerkin  (EBG)  methods,  such  as  FE,  SE, 
FV,  and  DG,  incorporate  the  orography  via  the  definition  of  the  grid.  EBG  methods  do  not  require  either 
orthogonal  grids  or  grids  with  specific  directions  (such  as  the  /  and  J  indices  in  FD  models);  EBG  models 
are  inherently  unstructured  and,  while  requiring  additional  data  structures  for  bookkeeping,  completely  liber¬ 
ate  the  method  from  the  grid.  This  freedom  from  the  grid  has  major  repercussions  in  the  implementation  of 
these  methods  on  distributed-memory  computers  in  that  no  halo  is  required  which  translates  into  truly  local 
algorithms  that  require  very  little  communication  across  processors;  instead,  the  communication  stencil  con¬ 
sists  of  the  perimeter  values  of  each  processor  (see  [21]).  The  advantage  that  SE  and  DG  methods  have  over 
the  FD  and  FV  methods  is  that  high-order  solutions  (greater  than  fourth  order)  can  be  constructed  quite  nat¬ 
urally  within  the  framework  -  such  high-order  properties  are  desirable  because  they  reduce  the  dispersion 
errors  associated  with  the  discrete  spatial  operators  [19].  In  fact,  the  SE  and  DG  formulations  proposed  in 
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this  paper  allow  for  arbitrarily  high-order  spatial  operators  to  be  constructed  by  an  input  parameter;  all  the 
results  presented  in  Section  6  use  either  8th  or  10th  order  polynomials  per  element. 

It  should  be  mentioned  that  there  are  other  numerical  methods  in  the  literature  that  have  much  to  offer 
nonhydrostatic  mesoscale  atmospheric  modeling.  Most  notably  are  the  weighted  essentially  non-oscillatory 
(WENO)  [15,27,45],  spectral  finite  volume  (SFV)  [56],  and  spectral  difference  (SD)  [36]  methods.  These  meth¬ 
ods  possess  many  of  the  best  features  of  the  DG  method  (such  as  high-order  accuracy,  conservation,  and  the 
promise  of  monotonicity).  The  reason  we  have  chosen  SE  and  DG  methods  for  our  study  is  that  they  can  be 
constructed  to  high-order  on  unstructured  grids  (either  with  quadrilateral  or  triangular  elements).  WENO, 
SFV,  and  SD  methods  do  not  offer  the  same  level  of  geometric  flexibility  to  use  unstructured  grids  in  combi¬ 
nation  with  high-order  accuracy.  There  are  triangular  WENO  schemes  but  only  for  very  specific  stencils  (see 
[27]  for  unstructured  WENO  methods  up  to  fourth  order);  this  is  true  also  for  the  SD  method  although  this 
method  is  still  in  its  infancy  and  higher  order  stencils  will  undoubtedly  be  constructed  in  the  future.  The  SFV 
method,  unfortunately,  requires  the  elements  to  be  planar  which  is  of  little  use  in  our  future  work.  In  this 
paper  we  only  consider  x-z  models  but  for  the  three-dimensional  model,  we  envision  using  curved  triangular 
prisms  in  order  to  take  advantage  of  unstructured  grids  in  the  horizontal  (x-y)  direction;  this  will  be  most  eas¬ 
ily  achieved  with  either  SE  or  DG  methods. 

The  remainder  of  the  paper  is  organized  as  follows.  Section  2  describes  the  three  forms  of  the  Euler  and 
Navier-Stokes  equations  that  we  study  in  this  work.  In  Section  3  we  discuss  the  seven  test  cases  used  to  com¬ 
pare  our  numerical  models.  Section  4  describes  the  spectral  element  and  discontinuous  Galerkin  formulation 
of  the  governing  equations  including  the  basis  functions,  numerical  fluxes,  and  boundary  conditions.  In  Sec¬ 
tion  5  we  describe  the  explicit  third  order  Runge-Kutta  method  we  use  to  march  the  equations  in  time  and  the 
filters  for  maintaining  stability.  In  Section  6  we  present  the  results  of  the  SE  and  DG  models  using  all  seven 
test  cases.  Finally,  in  Section  7  we  summarize  the  key  findings  of  this  research  and  propose  future  directions. 

2.  Governing  equations 

In  this  paper  we  study  three  different  forms  of  the  equations  that  govern  the  dynamics  of  nonhydrostatic  atmo¬ 
spheric  processes.  These  three  equation  sets  are  the  Euler  equations  that  have  been  used  for  many  years  in  com¬ 
putational  fluid  dynamics  (CFD).  One  of  the  sets  that  we  explore  is  the  complete  compressible  Navier-Stokes 
equations  with  the  physical  viscosity  defined  by  the  Stokes  hypothesis.  It  should  be  pointed  out,  however,  that 
we  use  viscosity  for  the  Navier-Stokes  equations  only  for  comparing  with  other  previously  published  model 
results  (see  the  density  current  test  in  Section  6);  the  main  focus  of  this  work  is  on  the  inviscid  form  of  the  equa¬ 
tions  (i.e.  the  Euler  equations).  Specifically,  we  study  the  following  three  equation  sets: 

1.  (set  1)  the  non-conservative  form  using  Exner  pressure,  momentum,  and  potential  temperature, 

2.  (set  2)  the  conservative  form  using  density,  momentum,  and  potential  temperature,  and 

3.  (set  3)  the  conservative  form  using  density,  momentum,  and  total  energy. 

For  the  purposes  of  this  study  we  restrict  ourselves  to  two  dimensions  (x-z)  and  omit  the  Coriolis  terms. 
2.1.  Equation  set  1 

Equation  set  1 ,  which  has  been  used  extensively  in  mesoscale  modeling,  reads 

071  „  R  „ 

—  +  U  ■  V 71  H - 7rV  •  u  =  0, 

at  cv 

du  t  . . , 

— -+  h  •  Vu  +  cp6\k  =  —gk  +  /A7“m,  (1) 

^+h-  \e  =  fV26, 

j-  /  \  r 

where  the  solution  vector  is  (7 i,*r  ,0)  ,  n  =  is  the  Exner  pressure,  u  =  (u,w)  is  the  velocity  field, 

0  =  ^  is  the  potential  temperature,  and  T  denotes  the  transpose  operator.  In  these  equations  P  is  the  pressure, 


3852 


F.X.  Giraldo,  M.  RestelH  I  Journal  of  Computational  Physics  227  (2008)  3849-3877 


P0  is  the  pressure  at  the  surface  (P0  =  1  x  105  Pa)  and  T  is  the  temperature.  Other  variables  and  symbols 
requiring  definition  are  the  gradient  operator  V  =  (T,T)T,  the  gravitational  constant  g,  the  gas  constant 
R  =  cp  —  cv,  the  specific  heats  for  constant  pressure  and  volume,  cp  and  c,„  the  dynamic  viscosity  //.  and  the 
directional  vector  along  the  vertical  (z)  direction  k  =  (0, 1)T.  The  viscosity,  p,  is  zero  for  all  tests  except  for 
the  density  current. 

The  advantage  that  set  1  has  over  the  other  two  sets  is  that  it  is  completely  self-contained,  that  is,  it  can  be 
solved  with  only  four  equations  defining  the  four  unknowns  (five  in  three-dimensions).  The  only  disadvantage 
is  that  a  model  based  on  these  equations  cannot  conserve  either  mass  or  energy.  Note  that  the  mass  equation  is 
defined  by  a  conservation-like  law  for  the  Exner  pressure  which  cannot  be  formally  conserved.  Existing  meso- 
scale  models  based  on  similar  equations  to  set  1  include  (but  are  not  limited  to)  ARPS  [60]  (University  of 
Oklahoma),  COAMPS  [26]  (US  Navy),  LM  [17]  (German  Weather  Service),  MM5  [25]  (Penn  State/NCAR), 
and  NMM  [29]  (NCEP). 

Introducing  the  following  splitting  of  the  Exner  pressure  and  potential  temperature  n(x,  t)  =  n{z)  +  7t'(jt,  t) 
and  0(x,t)  =  6(z)  +  &{x,t)  where  the  mean  values  are  in  hydrostatic  balance: 


-An 

c,8-=-g 


(2) 


allows  us  to  rewrite  Eq.  (1)  as 

An'  ,  An  R  ,  , 

— — b  u  ■  \n  +  w— — I - (n  +  n)\  ■  u  =  0. 

at  dr  cv 

0  q/ 

—  +  u  ■  \u  +  cp9Yn'  =  +  /<V2h, 

+  u  ■  \e'  +  w^-  =  pv2d, 

0?  dz 


(3) 


which  has  been  expanded  and  simplified  in  order  to  enforce  hydrostasis;  Eq.  (3)  is  the  form  used  for  all  the  test 
cases  in  Section  6. 


2.2.  Equation  set  2 


Equation  set  2  is  gaining  popularity  in  the  literature  because  it  is  not  too  dissimilar  from  set  1  and  is  in 
conservation  form  (for  the  inviscid  case  only).  These  equations  are  written  as  follows: 

|  +  V(p„)=0, 

V- (p«(g)«  +  PX2)  =  -pgk+  V  •  (np\u),  (4) 

^-+V-(p0u)  =  \-(pp\6), 

where  the  conserved  variables  are  (p,  puT ,  p0)T ,  p  is  the  density,  u  =  (m,w)t  is  the  velocity  field,  and  0  is  the 
potential  temperature  which  we  have  defined  previously.  The  pressure  P  which  appears  in  the  momentum 
equation  is  obtained  by  the  equation  of  state 


P  = 


(5) 


and  is  required  in  order  to  close  the  system.  Additional  terms  requiring  definition  are  the  tensor  product  ®  and 
the  rank-2  identity  matrix  X2;  this  term  essentially  converts  the  pressure  (which  is  a  scalar)  into  a  tensor. 

The  advantage  that  set  2  has  over  set  1  is  that  it  is  in  conservation  form,  which  when  using  methods  that  are 
formally  conservative,  allows  the  model  to  conserve  all  quantities.  Note,  however,  that  if  the  discretization 
method  is  not  formally  conservative,  then  this  set  should  have  little  or  no  advantage  over  set  1 .  Existing  meso- 
scale  models  based  on  this  equation  set  includes  WRF  [48]  and  the  model  by  Ahmad  and  Lindeman  [1], 
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Upon  splitting  of  the  density  and  potential  temperature  as  p(x,t)  =  p(z)  +  p'(x,t)  and 
d(x,  t )  =  9(z)  +  6'(x,  t)  where  the  mean  values  are  in  hydrostatic  balance,  Eq.  (4)  can  be  written  as 

^-+v.(P«)  =  0, 

+  V  •  (pu  ®  u  +  P': Zb)  =  -p'gk  +  V  •  (pp\u) ,  (6) 

d^-+\-(p0u)=\-(pp\6), 

where  P'  =  P  -  P  and  P  =  P(p,  9);  Eq.  (6)  is  the  form  used  for  all  the  test  cases  in  Section  6. 

2.3.  Equation  set  3 


Equation  set  3  is  the  form  used  in  computational  fluid  dynamics  (CFD,  e.g.  [18])  and  has  not  been  used  in 
atmospheric  studies  because  the  energy  equation  uses  total  energy  rather  than  potential  temperature  which 
then  requires  an  additional  step  to  compute  potential  temperature  in  order  to  use  existing  (moist)  sub-grid 
scale  physical  parameterizations.  However,  as  we  show  in  this  paper,  this  equation  set  has  some  advantages 
that  may  be  worth  considering  for  the  development  of  future  mesoscale  numerical  weather  prediction 
(NWP)  models. 

These  equations  are  written  as  follows: 

|+v-w  =  o, 

^  +  V  •  {pu  <g>  u  +  PI2)  =  -pgk  +  V  •  (7) 

^  +  V.[(pe  +  P)u]  =  \-Fr, 

where  the  conserved  variables  are  (p,puT,pe)T,  p  is  the  density,  u=(u,w)T  is  the  velocity  field, 
e  =  c,  T  +  [«■«  +  <p  is  the  total  energy,  and  cp  =  gz  is  the  geopotential  height.  The  pressure  P  is  obtained 
by  the  equation  of  state  which,  in  terms  of  the  solution  variables,  is  written  as 


P  = 


-uu-tp 


(8) 


Note  that  the  pressure,  Eq.  (8),  for  set  3  is  less  expensive  to  compute  than  the  pressure  for  set  2  (Eq.  (5)).  This 
will  be  shown  to  have  repercussions  in  the  relative  computational  costs  of  these  two  equation  sets. 

The  viscous  fluxes  FV1SC  are  defined  as  follows: 


FT  =  /i[V«  +  (v«f  +  2(V  •  «)Z2]  (9) 

and 

Z™c  =  u  ■  +  —\T,  (10) 

Pr 

where  y  =  cf  is  the  specific  heat  ratio,  X  =  —  |  comes  from  the  Stokes  hypothesis,  and  Pr  is  the  Prandtl  number. 

This  equation  set  directly  represents  the  fundamental  principles  of  conservation  of  mass,  momentum  and 
energy,  and  is  physically  consistent  both  in  the  inviscid  and  in  the  viscous  regimes.  In  particular,  when  viscos¬ 
ity  is  present,  it  naturally  accounts  for  the  dissipative  conversion  of  potential  and  kinetic  energy  into  internal 
energy.  In  addition,  from  a  more  computationally  oriented  viewpoint,  it  should  be  noted  that,  since  this  is  the 
equation  set  traditionally  employed  in  CFD,  much  of  the  machinery  developed  in  this  field  can  be  recycled; 
examples  include  total  variation  diminishing  (TVD)  schemes,  slope  limiters,  and  approximate  Riemann  solv¬ 
ers  (e.g.  see  [18,55]). 

Introducing  the  following  splitting  of  the  density  and  total  energy  p(x1t)  =  p(z)  +  p'(x,t)  and 
e(x,t)  =  e(z)  +e'(x,t),  where  the  mean  values  are  in  hydrostatic  balance,  allows  us  to  rewrite  Eq.  (7)  as 
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dp' 

i  +  v.(p„)  =  o, 

^  +  V  •  (pu  ®  «  +  P'1,)  =  -p'gk  +  V  •  F™c,  (11) 

^  +  V-[(pe  +  JP)«]  =  V-Fri 

which  is  the  form  used  in  all  the  test  cases  in  Section  6;  note  that  we  were  not  able  to  find  a  single  mesoscale 
model  in  the  literature  which  uses  this  equation  set. 

2.4.  Viscous  stresses 

It  should  be  mentioned  that  only  equation  set  3  has  the  true  viscous  stresses;  we  use  ad  hoc  diffusion  terms 
for  sets  1  and  2  merely  for  convenience.  The  viscous  terms  of  set  2  would  be  identical  to  those  of  set  3  for  the 
momentum  if 

v  •  (pp\u)  -  v  •  r;sc, 

where  the  arrow  signifies  that  the  term  on  the  left  should  be  replaced  by  the  one  on  the  right.  For  the  potential 
temperature,  after  much  algebra,  one  would  arrive  at 

These  viscous  terms  can  then  be  used  to  compute  the  true  stresses  of  set  1 .  Specifically,  for  set  1  we  would  get 
for  momentum 

/rV2«  ^ -U  •  F™c 

and  for  potential  temperature 

p\2e  ->  - — (\  ■  +  fvisc  •  Vu) . 

cpnp  \  \  Pr  J  ) 

There  are  two  reasons  why  we  do  not  use  these  stresses:  the  first  reason  is  that  for  the  purposes  of  intercom¬ 
parison  of  our  models  with  other  previously  published  models  we  adhere  to  the  equations  used  in  Straka  et  al. 
[51]  which  use  viscous  stresses  similar  to  those  that  we  define  for  sets  1  (Eq.  (1))  and  2  (Eq.  (4)).  The  second 
reason  for  not  using  these  stresses  (at  least  for  set  2)  is  that  the  equations  would  no  longer  be  in  conservation 
form.  Due  to  these  issues  alone  we  recommend  set  3  to  be  used  for  all  future  nonhydrostatic  atmospheric  mod¬ 
els.  For  the  purposes  of  this  study  it  is  not  so  important  to  use  the  true  viscous  stresses  for  sets  1  and  2  since 
viscosity  is  only  used  for  one  test  case.  Since  no  previous  work  using  set  3  was  found,  we  use  the  true  Navier- 
Stokes  viscosity  (Eqs.  (9)  and  (10)).  Due  to  this  discrepancy  in  the  viscous  stresses  between  the  three  equation 
sets,  we  expect  slight  differences  in  the  results;  the  question  is  how  large  will  these  differences  be. 

3.  Definition  of  the  test  cases 


We  now  describe  the  seven  test  cases  used  to  validate  the  numerical  models.  The  set  of  test  cases  that  we  use 
is  based  on  a  list  first  proposed  by  Skamarock  et  al.  [47];  we  have  added  two  bubble  tests  to  the  set.  None  of 
the  cases  uses  Coriolis  and  only  the  density  current  test  requires  viscosity. 

For  all  cases,  we  shall  define  the  initial  conditions  in  terms  of  Exner  pressure,  n,  and  potential  temperature, 
6.  However,  for  equation  sets  2  and  3  we  require  density,  p ,  and  for  set  3  total  energy,  e.  Thus  the  conversion 
from  7i  and  9  to  p  and  e  is  as  follows: 


and 
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e  =  cv07i  +  -u  ■  u  +  gz. 

To  simplify  the  description  of  the  initial  conditions  we  shall  write  the  Exner  pressure  and  potential  temper¬ 
ature  in  terms  of  a  mean  value  plus  a  perturbation,  that  is  n(x,  t)  =  n{z)  +  7t'(jc,  t)  and  0(x,  t)  =  0(z)  +  &(x,  t) 
where  the  reference  state  is  in  hydrostatic  balance. 


3.1.  Case  1:  Inertia- gravity  waves 


The  nonhydrostatic  inertia-gravity  wave  test  involves  the  evolution  of  a  potential  temperature  perturbation 
in  a  channel  with  periodic  boundary  conditions  on  the  left  and  right  boundaries.  The  initial  perturbation  radi¬ 
ates  to  the  left  and  right  symmetrically,  but  because  of  the  mean  horizontal  flow,  does  not  remain  centered 
about  the  initial  position.  The  initial  conditions  we  use  are  identical  to  those  of  Skamarock  and  Klemp 
[46];  we  provide  the  test  case  definition  below  for  completeness.  The  initial  state  of  the  atmosphere  is  taken 
to  have  a  constant  mean  flow  of  u  =  20  m/s  in  a  uniformly  stratified  atmosphere  with  a  Brunt-Vaisala  fre¬ 
quency  of  AT  =  0.01/s.  Using  the  definition  of  Brunt-Vaisala  frequency 

fC2  =  g±(\nd)  (12) 

yields 

e  =  60&,  (13) 

where  the  constant  of  integration  is  60  =  300  K.  The  Exner  pressure  is  obtained  from  the  hydrostatic  balance, 
Eq.  (2),  as 

S=l+— ^(e-^-l).  (14) 


In  addition,  the  potential  temperature  is  perturbed  by  the  amount 


O'  =  0C 


2  ! 


where  6C  =  0.01  °C,  hc  =  10,000  m,  ac  =  5000  m,  xc  =  100,000  m,  and  nc  =  3.14159265  is  the  Archimedes’ 
(trigonometric)  constant. 

The  domain  is  defined  as  (x.z)  e  [0, 300, 000]  x  [0, 10, 000]  m  with  t  e  [0,  3000]  s.  No-flux  boundary  condi¬ 
tions  are  used  along  the  bottom  and  top  boundaries  while  periodic  boundary  conditions  are  used  along  the 
lateral  boundaries. 


3.2.  Case  2:  Rising  thermal  bubble 


The  rising  thermal  bubble  test  shows  the  evolution  of  a  warm  bubble  in  a  constant  potential  temperature 
environment.  Because  the  bubble  is  warmer  than  the  ambient  air  it  rises  while  deforming  as  a  consequence  of 
the  shearing  motion  caused  by  the  velocity  field  gradients  until  it  forms  a  mushroom  cloud.  The  initial  con¬ 
ditions  we  use  are  similar  to  those  of  Robert  [42]  with  the  air  mass  initially  at  rest  and  in  hydrostatic  balance. 
Setting  0  =  constant  in  (2)  immediately  provides  for  the  Exner  pressure.  To  drive  the  motion  of  the  air,  the 
following  potential  temperature  perturbation  is  then  added: 


f° 

1  +COS  (f)' 

for  r  >  rc. 
for  r  if  rc, 


where  9C  =  0.5  °C,  nc  is  the  trigonometric  constant,  r  =  \J (x-  xc)2  +  (z  —  zc )2  with  the  following  constants: 
0  =  300  K,  rc  =  250  m,  and  (x,z)  e  [0, 1000]2  m  with  t  €  [0, 700]  s  and  (xc,zc)  =  (500,  350)  m.  The  boundary 
conditions  for  all  four  boundaries  are  no-flux. 
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3.3.  Case  3:  Robert  smooth  bubble 


Test  case  3  is  similar  to  the  Robert  smooth  bubble  presented  in  [42];  however,  we  use  a  cosine  wave  instead 
of  the  Gaussian  used  in  [42],  The  definition  of  this  case  is  almost  identical  to  the  rising  thermal  bubble  except 
that  the  domain  is  now  defined  as  (x,z)  £  [0, 1000]  x  [0, 1500]  m  with  t  £  [0, 800]  s  and  (xc,zc)  =  (500, 260)  m. 
Recall  that  in  case  2,  the  bubble  is  initially  positioned  at  yc  =  350  whereas  now  we  lower  this  position  and  let 
the  bubble  evolve  for  an  additional  100  s  (from  700  s  in  case  2  to  800  s  in  case  3).  The  longer  integration 
requires  us  to  move  the  top  boundary  in  order  to  avoid  having  the  bubble  hit  this  boundary. 


3.4.  Case  4:  Density  current 


The  density  current  test  was  proposed  in  [51]  and  concerns  the  evolution  of  a  cold  bubble  dropped  in  a  neu¬ 
trally  stratified  atmosphere.  Because  the  bubble  is  cold,  it  sinks  eventually  hitting  the  ground.  At  this  point,  the 
bubble  begins  to  shear  as  it  travels  along  the  ground  forming  Kelvin-Helmholtz  rotors.  As  discussed  in  [51], 
viscosity  is  required  in  order  to  obtain  a  grid-converged  solution.  The  initial  conditions  for  this  case  are  quite 
similar  to  those  of  the  rising  thermal  bubble;  however,  the  differences  are  in  the  domain  size  and  the  shape  of 
the  cold  cosine  bubble.  The  potential  temperature  perturbation  is  defined  as 

0'  =  y  [1  +  cos(7tcr)], 
where  6C  =  — 15  °C 


and  rc=  1.  The  domain  is  defined  as  (x,z)  £  [0, 25, 600]  x  [0, 6400]  m  with  t  £  [0, 900]  s  and  the  center  of  the 
bubble  is  originally  positioned  at  (xc,zc)  =  (0, 3000)  m  with  the  size  of  the  bubble  defined  by  the  parameters 
(xr,zr)  =  (4000, 2000)  m.  The  boundary  conditions  for  all  four  boundaries  are  no-flux  and  a  dynamic  viscosity 
of  ju  =  75  m2/s  is  used.  Note  that  we  have  defined  only  half  the  domain  in  the  horizontal  as  proposed  in  Stra- 
ka  et  al.  [51], 


3.5.  Case  5:  Schar  mountain 


The  Schar  mountain  test  concerns  the  steady-state  solution  of  hydrostatic  flow  over  a  five-peak  mountain 
chain,  with  steady  inflow  and  outflow  boundary  conditions.  The  initial  conditions  and  mountain  profile  are 
given  in  Schar  et  al.  [43]  which  we  now  provide  for  completeness.  The  initial  state  of  the  atmosphere  is  taken 
to  have  a  constant  mean  flow  of  u  =  10  m/s  in  a  uniformly  stratified  atmosphere  with  a  Brunt-Vaisala  fre¬ 
quency  of  Af  =  0.01/s.  Using  the  definition  of  Brunt-Vaisala  frequency,  Eq.  ( 12)  yields  the  reference  potential 
temperature  given  in  Eq.  (13),  the  Exner  pressure  given  in  Eq.  (14)  where  the  constant  of  integration  is  taken 
to  be  do  =  280  K.  The  domain  is  defined  as  (x,z)  €  [—25,000,25,000]  x  [0,21,000]  m  with  t  £  [0, 10]  h.  The 
mountain  is  defined  as 


with  the  parameters  hc  =  250  m,  Xc  =  4000  m,  and  ac  =  5000  m.  This  profile  is  shown  in  Fig.  1  where  the  axes 
have  been  magnified  as  follows  (x,z)  £  [—10, 000, 10, 000]  x  [0, 1000]  m  in  order  to  show  the  five  peaks  of  the 
mountain.  No-flux  boundary  conditions  are  used  along  the  bottom  surface  while  non-reflecting  boundary  con¬ 
ditions  are  used  along  the  top  and  lateral  boundaries.  Note  that  this  test  is  in  the  hydrostatic  range  since  ^fL>  1  • 


3.6.  Case  6:  Linear  hydrostatic  mountain 


The  linear  hydrostatic  mountain  test  involves  the  steady-state  solution  of  linear  hydrostatic  flow  over  a 
single-peaked  mountain  with  constant  inflow  and  outflow  boundary  conditions.  The  initial  conditions  and 
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Fig.  1.  Case  5:  The  Schar  Mountain  profile  magnified  such  that  (x,z)  £  [—10,000, 10,000]  x  [0, 1000]  m. 


mountain  profile  are  presented  in  Smith  [49]  and  Durran  and  Klemp  [11],  The  initial  state  of  the  atmosphere 
consists  of  a  constant  mean  flow  with  w  =  20  m/s  in  an  isothermal  atmosphere  with  constant  temperature 
T  =  250  K.  Since  in  the  isothermal  case  we  have  A/"  =  — £==,  the  potential  temperature  and  Exner  pressure  pro¬ 


files  are  immediately  provided  by  Eqs.  (13)  and  (14),  this  latter  reducing  to 


7t  =  e  cpt  . 

The  domain  is  defined  as  follows:  (x,z)  e  [0,240,000]  x  [0,30,000]  m  with  t  g  [0, 10]  h  with  the  versiera  di 
Agnesi  mountain  profile  defined  as 

h(x,z)  = - hc  (15) 

1  +  (t?) 

where  hc  =  1  m,  xc  =  120, 000  m,  and  ac  =  10, 000  m.  It  can  be  verified  that  J^fL  >  1,  so  that  the  flow  is  in  the 
hydrostatic  range.  Non-reflecting  boundary  conditions  are  required  along  the  lateral  and  top  boundaries.  In 
addition,  no-flux  boundary  conditions  are  used  along  the  bottom. 


3. 7.  Case  7:  Linear  nonhydrostatic  mountain 

The  linear  nonhydrostatic  mountain  test  involves  the  steady-state  solution  of  linear  nonhydrostatic  flow 
over  a  single-peaked  mountain  with  constant  inflow  and  outflow  boundary  conditions.  For  this  test  the  initial 
state  of  the  atmosphere  consists  of  a  constant  mean  flow  of  u  =  10  m/s  in  a  uniformly  stratified  atmosphere 
with  a  Brunt-Vaisala  frequency  of  J\f  =  0.01  /s.  Potential  temperature  and  Exner  pressure  profiles  are  given  by 
Eqs.  (13)  and  (14),  respectively,  with  00  =  280  K.  The  mountain  profile  given  by  Eq.  (15)  is  used,  with  the 
numerical  parameters  now  specified  as  hc  =  1  m,  xc  =  72,000  m,  and  a,  =  1000  m.  The  computational 
domain  is  (x,  z)  g  [0,144,000]  x  [0,30,000]  m  with  tg  [0,5]  h.  Note  that  since  AA  =  1  the  flow  is  then  in 
the  nonhydrostatic  regime.  No-flux  boundary  conditions  are  used  along  the  bottom  and  non-reflecting  bound¬ 
ary  conditions  along  the  lateral  and  top  boundaries. 


4.  Spatial  discretization  of  the  governing  equations 

For  the  sake  of  brevity,  we  shall  only  discuss  the  discretization  of  set  3  which  we  now  write  in  the  following 
compact  vector  form 
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dq 

dt 


\-F  =  S(q), 


(16) 


where  q  =  (p,  UT ,E)r  is  the  solution  vector  with  U  =  pu  and  E  =  pe,  F  =  -  FV1SC  is  the  total  flux  where 

l  U  \ 


F™(q)=  ^  +  PT 2 

(E+P)U 
p 

is  the  inviscid  flux  tensor, 


V 


(17) 


0 


F™c(q)  =  /(  |  V«  +  (V«)  +  1(V  •  li)I2 

u  [\u  +  (Vuf  +  ■  u)I2]  +%\T / 

is  the  viscous  flux  tensor,  and 


(18) 


(19) 


S(q)  =  - 

V  o  J 

is  the  source  function.  Eq.  (16)  must  be  solved  in  Q  x  (0,  7>inal),  with  Q  c  IR2. 

4.1.  Basis  functions,  differentiation,  and  integration 

To  define  the  local  operators  which  shall  be  used  to  construct  the  global  approximation  of  the  solution  we 
begin  by  decomposing  the  domain  Q  into  Ne  non-overlapping  quadrilateral  elements  such  that 

Ne 

Q=\jQe- 

e—\ 

We  now  define  the  reference  element  I  =  [ —  1 ,  l]2  and  introduce  for  each  element  Q,.  the  smooth,  bijective  trans¬ 
formation  Tqb  such  that  Qe  =  ^^(l).  The  notation  x  =  (Faf  f)  will  also  be  used,  with  x  =  (x,z)  €  Qe  and 
§  =  G  I .  Associated  with  the  map  TQe  is  the  Jacobian  Jor  =  with  determinant  Jq„  .  The  Jacobian  deter¬ 
minant  of  To,,  restricted  to  the  boundary  of  I  is  denoted  by  JSQ  .  In  both  the  SE  and  DG  methods,  an  approxima¬ 
tion  qN  of  the  solution  of  Eq.  (16)  is  sought  such  that,  at  each  time  level  t  e  [0,  rfinai] 

e  =  1, . . .  ,Ne, 

where  if  is  a  function  of  the  space  VN(\)  of  bivariate  polynomials  of  degree  lower  than  or  equal  to  N  in  i 
VN{\)  =  Span^VI  m,n  <  A,  (£,r,)  el},  1. 

A  critical  step  is  now  the  definition  of  a  basis  {^k}f=l,E  =  (N  +  l)2,  for  VN{\).  The  tensor  product-structure  of 
I  allows  us  to  construct  such  a  basis  as 

<M$)  =  h‘^)hM), 

where  {hi}f=0  is  a  basis  for  VN{[- 1, 1])  and  an  index  kij  is  biunivocally  associated  with  the  pair  (i,j),  with 
ij  =  0, . . . ,  N.  The  basis  functions  hff)  are  in  turn  constructed  as  the  Lagrangian  polynomials  associated  with 
the  Legendre-Gauss-Lobatto  (LGL)  points  c,  given  as  the  roots  of 

(i-aw=o, 

where  Pn(£)  is  the  Ath  order  Legendre  polynomial.  Notice  that  associated  with  these  points  are  the  Gaussian 
quadrature  weights 


a>i  = 


A  (A  +  1)  VVfe) 


1 


2 
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This  fact  will  be  exploited  momentarily  for  the  evaluation  of  local  element  integrals  of  qN. 

With  these  definitions,  we  have 

K 

t)\Qe  =  'l'k{Fa](x))qk(t),  e=l,-..  ,Ne,  (20) 

k=  1 


where  we  introduce  the  grid  points  xkij  =  ((£,■,  if))  and  the  grid  point  values  qk(t)  =  qN(  xk,  t).  The  expan¬ 

sion  given  in  Eq.  (20)  is  essential  in  order  to  compute  local  element-wise  derivatives  and  integrals.  Concerning 
derivatives,  in  fact,  it  immediately  provides 


q9n 

dx 


(x,t) 


=  E  4z[U^(x))]qt(t),  ^(*,0  =  w- 


tie  k=  1  k=  1 

Concerning  the  computations  of  integrals,  the  expansion  defined  by  Eq.  (20)  yields 

[  qN(x,t)dx  =  I  qN(x{!;),t)JQe{Z)di;~^(OiCOjqklj{t)\Jae(Zi,rij)\. 

J  Qe  j  I  iJ—0 


(21) 


(22) 


4.2.  Spectral  element  method 

In  the  spectral  element  method,  the  following  weak  form  approximation  of  Eq.  (16)  is  considered:  find 
qN{-,  t)  G  VSE  such  that 

f{9n) ) dQ= L  v</>  e  f"e’  (23) 

where 

Vf  =  {</>  G  H\Q)  :  (t>\Qe  =  .A  o^1,  with  ,//  G  7V(0,e=  1,. . .  ,2Ve}.  (24) 

Notice  that  the  two  requirements  <j)  G  Hl(Q)  and  (f>\Qe  =  1/70  Tq]  imply  VSE  C  C°(£2).  Integrating  Eq.  (23)  by 
parts  yields  the  equivalent  problem:  find  qN(-,t)  G  Vff  such  that 

J^dQ  +  J^n-  F(qN)dr  -  j  V</>  •  .F(-7w)df2  =  J  <l>S(qN)dQ,  G  Vf,  (25) 

where  //  is  the  outward  pointing  normal  vector  on  the  domain  boundary  T.  We  point  out  that,  although  Eqs. 
(23)  and  (25)  are  in  principle  equivalent,  they  yield  different  numerical  formulations  when  approximate  quad¬ 
rature  formulas  are  adopted  for  the  computation  of  the  flux  integrals,  as  is  often  the  case.  Such  quadrature 
formulas  in  fact  increase  the  computational  efficiency  of  the  scheme  and  diagonalize  the  mass  matrix  M  to 
be  defined  shortly.  In  this  paper,  we  shall  use  the  form  (25),  since  it  allows  for  global  conservation  of  mass, 
momentum,  and  energy  even  when  approximating  the  integrals  as  in  (22)  [52], 

By  virtue  of  Eqs.  (21)  and  (22),  Eq.  (25)  yields  the  matrix  problem 

J  +  ( M°)TF(q )  -  DTF(q)  =  S(q),  (26) 

where 

MS=M~XMS,  D  =  MlD 

and  M,  Ms  and  D  are  the  global  mass,  boundary  and  differentiation  matrices,  respectively.  These  latter  are  in 
turn  constructed  from  their  local  counterparts  Me,Ms,e  and  De  by  means  of  the  direct  stiffness  summation 
procedure 

Ne  Ne  Ne 

M  =  /\Me ,  Ms  =  f\Ms'e,  D=/\De , 

e=l  e=l  e=l 
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where  denotes  mapping  the  local  degrees  of  freedom  of  each  element  Qe  to  the  corresponding  global  de¬ 
grees  of  freedom  in  Q  followed  by  summation  (see  [21]  for  further  details).  The  local  matrices  are  then  defined 
as 

Mhk  =  ™Vne(Sh)\SHk,  Dlk  =  w*I-/q.(«*)|V^(x*),  M\l  =  ^h\JsQe^h)\dhkn{xh): 

where  h,k  =  1 ,K,dhk  is  the  Kronecker  delta,  ty) ,  wkij  =  cOjCOj ,  w\r  =  m,  for  j  =  0  or  j  =  N  and 

vfk  =  coj  for  /  =  0  or  i  =  N. 

4.3.  Discontinuous  Galerkin  method 


The  discontinuous  Galerkin  method  comes  in  various  forms  and  flavors.  For  example,  one  could  use  either 
modal  (amplitude-frequency  space)  or  nodal  (physical  space)  polynomial  expansions.  In  this  paper,  we  exploit 
the  Lagrangian  basis  introduced  in  Section  4. 1  to  define  a  nodal  DG  formulation  with  inexact  integration  as  pro¬ 
posed  in  [20].  Thus,  the  DG  formulation  that  we  describe  below  is  essentially  the  discontinuous  version  of  the  SE 
formulation  illustrated  in  Section  4.2.  Alternatively,  one  could  use  exact  integration  (as  in  [23])  but  as  shown  in 
[20]  this  would  be  more  computationally  expensive  and  thereby  would  not  be  competitive  with  the  SE  method. 

In  the  DG  method,  Eq.  ( 16)  is  multiplied  by  a  test  function  <p  and  integrated  over  a  generic  element  Qe,  and 
the  exact  solution  is  replaced  by  its  approximation,  yielding 

L 0  (^r + v  ■ F^) dQ= ja  (27) 

where  qeN  denotes  the  degrees  of  freedom  collocated  in  Qe.  Applying  now  integration  by  parts  and  introducing 
the  numerical  flux  F* ,  the  following  problem  is  obtained:  find  qN(-,t )  €  F°G  such  that  VQ,..  e  =  1, . . . . 

[  MdQe+[  $n-F*(qN)dre-  [  V<j>  ■  F(<fN)dC2e  =  [  <j>S(qeN)dQ,  (28) 

JQe  Jre  JQe  JQe 

where 

k°G  =  {0  G  L\Q)  :  <£|Qc  =  tfr  o  with  +  G  VN(\),e=  1  ,...,Ne}.  (29) 

Notice  that,  contrary  to  Eq.  (24),  there  is  no  global  continuity  requirement,  so  that  F°G  ^  C°(Q).  This  is  pos¬ 
sible  because  in  Eq.  (25)  differentiability  is  required  separately  within  each  element,  and  not  within  the  entire 
domain  Q.  The  coupling  between  neighboring  elements  is  then  recovered  through  the  numerical  flux  F*,  which 
is  required  to  be  a  single  valued  function  on  the  interelement  boundaries  and  the  precise  definition  of  which  is 
given  below. 

By  virtue  of  Eqs.  (21)  and  (22),  Eq.  (28)  can  be  written  in  the  matrix  form 

^  +  (Ms’e)T F* (q)  -  {D'fFtf)  =  S(qe)  (30) 

where  =  (MeylMs-e  and  De  =  ( MeylDe .  Note  that  this  equation  can  be  simplified  to  yield  the  following 
semi-discrete  weak  form: 


dt 


(DlfF* 


■  SI  - 


Vs,' 


V\  Je, 


[mrF* 


(31) 


Further  integration  by  parts  in  Eq.  (28),  and  simplifying,  yields  the  following  semi-discrete  strong  DG 
formulation: 


dt 


-WuyF*+s*+ 


v\js,  i 

Wfl  Jei 


<ir-F*)t, 


(32) 


which  looks  very  similar  to  a  collocation  penalty  method  in  that  no  global  matrix  system  needs  to  be  con¬ 
structed;  it  should  be  understood  that  this  simple  form  is  only  possible  by  choosing  the  inexact  integration 
DG  form  introduced  in  [20]  which  uses  the  same  inexact  integration  theory  and  machinery  used  in  the  classical 
SE  method. 
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4.3.1.  Inviscid  numerical  flux 

Substituting  the  boundary  flux  (Eq.  (27))  with  the  numerical  flux  (Eq.  (28))  is  necessary  in  order  to  resolve 
the  ambiguity  in  the  definition  of  qN  on  1).  (where  it  can  attain  two  different  values)  and,  by  doing  this,  restor¬ 
ing  the  coupling  between  neighboring  elements.  Various  choices  are  possible  for  the  definition  of  F*;  in  the 
present  study  we  use  the  Rusanov  flux,  mainly  because  of  its  simplicity  but  also  because  we  have  detected 
no  discernible  differences  between  the  solution  obtained  with  this  flux  function  and  the  EILL,  EILLC,  or  Roe’s 
solver  (see  [18]).  It  is  quite  possible  that  for  more  stringent  test  cases  (e.g.  for  shocks),  differences  between  the 
various  Riemann  solvers  may  occur. 

Upon  arbitrarily  associating  to  each  interelement  boundary  an  intrinsic  tangential  direction,  so  that  a  left 
(L)  element  and  a  right  (R)  element  are  identified,  the  Rusanov  numerical  flux  can  be  expressed  as 

FtMv(q)  =\[FN(qLN)  +  FN(q R)  -  |A|(?R  -  </>],  (33) 

where  1  =  max(|UL|  +  y/a^,  |UR|  +  Va*)  with  UL’R  =  if'11  ■  n  being  the  normal  component  of  velocity  with  re¬ 
spect  to  the  edge  Te  and  a  =  y/yRT  being  the  speed  of  sound. 


4.3.2.  Viscous  numerical  flux 

For  the  viscous  form  of  the  equations,  the  DG  method  is  applied  in  a  slightly  different  way.  In  fact,  we  use 
the  local  discontinuous  Galerkin  method  first  proposed  in  [3]  and  later  analyzed  in  [8].  The  idea  is  very  similar 
to  ideas  put  forth  in  the  mixed  finite  element  approach.  In  other  words,  rather  than  discretizing  the  Laplacian 
directly  via  a  finite  element  approximation,  we  proceed  in  the  following  two-step  approach.  In  Eqs.  (9)  and 
(10)  we  define  the  viscous  fluxes  for  the  Navier-Stokes  equations  as 

U  =  \u  and  T  =  \T. 


Then,  we  apply  the  DG  discretization  in  weak  form  and  dividing  through  by  the  mass  matrix  obtaining 


u]  =  -TO* 


vo’vr 

Hf \j°\ 


«r«;  and  T*  =  -(D‘)rt 


wfl  JT 


(34) 


with  u*  =  I  (kl  +  mr)  and  T*  =  \  (TL  +  TR).  Note  that  other  choices  are  possible  but  we  have  chosen  for  sim¬ 
plicity  the  average  value  at  the  interface  (see  [8]  for  other  choices).  At  this  point,  we  can  now  discretize  the 
viscous  equations  as 


^  =  (be)T  Fe. 
d  t  1  ,J>  J 


%n- 


ns,’eF* 


where  Fe  =  F'm(qf  -  Fnsc(Ue,  Te),  F*  =  Fmv(q*)  -  Fvisc(U* ,T*),  U*  =\{UL 


(35) 

UR),  and  T*  =  I(TL  +  TR). 


4.4.  Boundary  conditions 


An  important  component  of  a  mesoscale  atmospheric  model  is  the  boundary  conditions.  For  the  test  cases 
considered  in  this  paper  (and  assuming  a  rectangular  domain)  the  boundary  conditions  are:  no-flux  along  the 
bottom  boundary  and  either  periodic,  no-flux,  or  non-reflecting  for  the  left  and  right  boundaries,  and  either 
no-flux  or  non-reflecting  for  the  top  boundary. 


4.4.1.  Non-reflecting  boundary  conditions 
Let  us  rewrite  the  governing  equations  as 


dq 
d  t 


- 


where  the  ellipsis  denotes  the  usual  terms  in  the  governing  equations  and  the  term  t  is  the  coefficient  of  the 
absorbing  sponge  layer  which  relaxes  the  numerical  solution  to  the  prescribed  reference  value  q.  In  terms 
of  an  explicit  time-integrator,  this  becomes  nothing  more  than 

qsponeed  =  q-<q-q), 
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where  q  is  the  value  obtained  from  the  dynamical  model,  q  is  the  reference  value,  and  q  spont'cd  is  the  sponged 
variable;  the  parameter  x  is  a  Lagrange  multiplier  that  is  found  using  a  similar  procedure  to  that  presented  in 
[11]. 

4.4.2.  No-flux  boundary  conditions 

The  no-flux  boundary  conditions  are  enforced  by  virtue  of  the  statement 

it  ■  u  =  0 

at  the  boundaries.  Thus,  we  seek  to  eliminate  the  normal  component  of  the  velocity  to  the  no-flux  boundary 
without  altering  the  tangential  component  (we  assume  free  slip  boundary  conditions  for  all  test  cases).  The 
tangent  vector  to  a  boundary  is  obtained  by  t  =  k  x  it  which  is  equal  to  t  =  —nzi  +  nxk.  Thus  we  solve  the  fol¬ 
lowing  2x2  system: 

/  nx  \  /  m  \  _  /  0  \ 

V  — «z  nxJ\wJ  \  uT  )  ’ 

where  uT  =  t  ■  u  is  the  tangential  component  of  velocity. 

For  the  inviscid  equations  (i.e.  the  Euler  equations)  the  no-flux  boundary  condition  on  the  momentum  is  all 
that  is  needed  for  all  three  equation  sets.  Concerning  the  viscous  case,  for  equation  sets  1  and  2,  we  simply 
impose  a  zero  viscous  flux  for  both  momentum  and  temperature.  For  equation  set  3,  however,  while  imposing 
a  zero  viscous  momentum  flux  is  still  appropriate,  enforcing  a  zero  viscous  energy  flux  would  result,  in  the  case 
of  a  neutral  atmosphere,  in  the  development  of  a  thermal  boundary  layer  at  the  top  and  bottom  boundaries.  In 
this  case  in  fact,  since  the  viscous  flux  is  controlled  by  the  temperature  gradient  rather  than  by  the  potential 
temperature  gradient,  the  viscous  flux 

pvisc  _  licp  di1  , 

e  “  Pr  dz 

is  required  in  order  to  maintain  the  neutral  stratification.  Hence,  this  latter  flux,  being  required  by  the  equa¬ 
tions  themselves  and  not  by  the  numerical  discretization,  will  be  imposed  in  both  the  SE  and  DG 
formulations. 

4.5.  Numerical  models 

In  Section  2  the  following  three  equation  sets  were  discussed:  set  1  given  by  Eq.  (3),  set  2  given  by  Eq.  (6), 
and  set  3  given  by  Eq.  (11).  The  eigenvalues  of  all  three  equation  sets  are  {U,  U ,  U  —  x/d.  U  +  fla)  where 
U  =  it  u  with  n  referring  to  a  unit  vector,  and  a  is  the  speed  of  sound.  It  is  important  to  know  the  eigenvalues 
in  order  to  choose  the  maximum  time-step  which  will  allow  the  model  to  run  in  a  stable  fashion,  and  for  con¬ 
structing  the  numerical  flux  in  the  discontinuous  Galerkin  models. 

In  Section  4  the  spectral  element  (SE)  and  discontinuous  Galerkin  (DG)  methods  were  described.  Based  on 
the  above  three  equation  sets  and  the  two  numerical  methods  we  have  developed  five  numerical  models  which 
we  compare  in  Section  6.  Table  1  summarizes  the  five  numerical  models  that  we  study. 

Table  1  shows  that  there  is  no  DG  model  for  equation  set  1;  this  is  because  equation  set  1  is  not  in  conser¬ 
vation  form.  Note  that  the  momentum  and  energy  equations  can  be  written  in  conservation  form  as  in  equa¬ 
tion  set  2.  However,  to  write  the  Exner  pressure  equation  in  conservation  form  requires  adding  (and 


Table  1 

Numerical  models  constructed  in  this  study  based  on  equation  sets  1.  2,  and  3  and  the  spectral  element  and  discontinuous  Galerkin 
methods 


Spectral  element 

Discontinuous  Galerkin 

Set  1  (Eq.  (3)) 

SE1 

- 

Set  2  (Eq.  (6)) 

SE2 

DG2 

Set  3  (Eq.  (11)) 

SE3 

DG3 
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subtracting)  terms  which  will  lead  to  a  source  term  that  cannot  be  avoided.  If  this  is  not  troubling  enough,  the 
fluxes  resulting  from  this  system  are  non-physical  and  hence  the  DG  method  will  not  work.  This  then  leaves  us 
with  only  five  numerical  models. 

5.  Time-integrator  and  filters 


In  order  to  advance  the  solution  in  time  while  retaining  some  level  of  high-order  accuracy,  we  use  a  strong 
stability  preserving  (SSP)  Runge-Kutta  third  order  (RK3)  method  of  the  type  first  proposed  by  Cockburn  and 
Shu  [9], 

Definition  5.1.  The  solution  q  is  said  to  be  strongly  stable  if  II?"IItv  ^  \W'+X IItv'^  n- 


Definition  5.2.  The  TV  norm  is  defined  as  (see,  e.g.  [54])  ||^r||xv  =  Y^!i=\  '  l</,+i  —  |. 


Let  us  write  the  semi-discrete  (in  space)  equations  as  follows 


dq 
d  t 


where  R(q)  is  the  SE  or  DG  representation  of  the  term  S(q)  —  V  •  F(q)  in  Eq.  (16)  (where  we  assume  a  method 
of  lines  approach).  The  SSP  temporal  discretization  of  this  vector  equation  is 

for  k  =  1, . . . ,  S  :  qk  =  ak0q"  +  akiqk~l  +  ft  A tR{qk~l), 

where  q°  =  q",qs  =  qn+l,  and  S  denotes  the  number  of  stages  used.  For  all  of  our  simulations  we  use  the 
SSPRK  (5, 3)  method  of  Spiteri  and  Ruuth  [50]  which  is  a  5-stage  third-order  strong  stability  preserving  Run¬ 
ge-Kutta  method. 

The  time-step  for  all  our  simulations  are  chosen  based  on  the  maximum  Courant  number  allowable  by  the 
time-integrator.  For  the  purposes  of  this  study  we  define  the  Courant  number  as 

fCAt\ e 

Courant  number  =  max  I  — —  )  Ve  €  [1 , . . . ,  Ne\ , 

\  As  J  ho 


where  C  =\U  +  \fti\  is  the  characteristic  speed,  U  =  n  u  is  the  velocity  in  the  direction  n,  a  is  the  sound  speed, 
and  A,v  =  V Ax2  +  Ar2  is  the  grid  spacing.  For  all  the  results  presented,  the  Courant  number  is  taken  to  be  1.3 
for  the  strong  form  SE  method,  Eq.  (23),  and  0.85  for  the  DG  method  (either  strong  or  weak),  and  for  the 
weak  form  SE  method,  Eq.  (25);  it  turns  out  that  imposing  the  boundary  conditions  weakly  reduces  the  max¬ 
imum  allowable  time-step.  We  shall  report  more  on  this  in  future  work. 

In  order  to  maintain  stability  for  long-time  integrations,  it  is  important  not  only  to  adhere  to  the  CFL  con¬ 
dition  but  also  to  control  the  Gibbs  oscillations  which  affect  all  high-order  methods  in  the  absence  of  a  dis¬ 
sipative  mechanism  (e.g.  viscosity,  slope  limiters,  etc.).  The  standard  approach  to  avoid  such  instabilities  is 
through  the  use  of  spatial  filters  of  the  Boyd-Vandeven  type  [6].  The  filtering  procedure  is  applied  after  each 
time-step  as  follows 

=  F9n- 

where  qFN  denotes  the  filtered  solution  of  the  state  variable  qN,  and  F  is  the  filter  matrix  defined  in  [21];  note  that 
the  filter  matrix  is  very  weak  -  in  fact,  only  the  highest  modes  are  reduced  by  5%. 


6.  Results 


In  this  section  we  validate  the  five  SE  and  DG  numerical  models  on  the  test  case  suite  of  seven  problems. 
For  four  of  the  seven  problems  that  do  not  have  analytic  solutions  we  use  symmetry  and  a  comparison  of 
extrema  as  metrics  to  discern  the  quality  of  the  models. 

For  the  three  mountain  problems  we  have  senri-analytic  solutions  based  on  linear  theory  which  require  the 
application  of  a  Fourier  transform  (see  [49]  for  details).  To  obtain  this  semi-analytic  solution  we  use  a  grid  of 
4000  (along  x)  by  100  (along  z)  points  in  order  to  ensure  that  the  solution  is  exact  within  machine  precision.  In 
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addition,  since  the  grid  for  the  numerical  and  analytic  solutions  are  different,  we  then  interpolate  both  solu¬ 
tions  to  a  grid  of  400  (in  x)  by  100  (in  z)  defining  a  specific  portion  of  the  domain.  This  portion  of  the  domain 
is  the  one  that  we  use  to  plot  and  compute  the  root-mean-square  (RMS)  errors.  We  define  the  root-mean- 
square  error  as 


\ 


^  ^  ^numerical  ^analytic  ^ 


z=l 


where  Np  =  40, 000. 

Following  Smith  [49],  for  the  hydrostatic  and  nonhydrostatic  mountains  we  define  the  momentum  flux  as 

/+oo 

p(z)u(x,  z)w(x,  z)  dx, 

•oo 

where  p(z)  is  the  reference  density  as  a  function  of  height.  From  linear  theory,  the  analytic  hydrostatic 
momentum  flux  is  given  as 

mH(z)  =  —  (36) 

where  the  superscript  H  signifies  Hydrostatic,  ps  and  us  are  the  reference  density  and  horizontal  velocity  values 
at  the  surface,  A f  is  the  Brunt-Vaisala  frequency,  and  hc  is  the  height  of  the  mountain. 

Following  Klemp  and  Durran  [32],  we  define  the  analytic  nonhydrostatic  momentum  flux  as 

wNH(z)  =  — 0.457/mh(z), 

where  mH(z)  is  given  by  Eq.  (36).  Note  that,  as  pointed  out  by  Klemp  and  Durran  [32],  this  solution  is  valid 
only  for  =  1,  which  we  specify  a  priori.  The  normalized  momentum  flux  that  we  present  below  is  then  de¬ 
fined  as  either  m(z)'mH(z)  or  m(z)'  wNH(z),  depending  on  whether  the  mountain  problem  is  in  the  hydrostatic 
or  nonhydrostatic  regime. 


6.1.  Case  1:  Inertia-gravity  waves 

Fig.  2a  shows  the  potential  temperature  perturbation  contours  for  DG3  after  3000  s  and  Fig.  2b  shows  the 
one-dimensional  profile  along  z  =  5000  nr  for  all  five  models.  Fig.  2b  shows  that  all  five  models  yield  identical 
solutions  -  this  is  especially  of  interest  since  the  models  use  different  numerical  methods  and  equation  sets.  The 
second  observation  is  that  the  profiles  are  perfectly  symmetric  about  the  position  x  =  160,000  m.  Note  that 


Fig.  2.  Case  1:  Inertia-gravity  waves.  Potential  temperature  perturbation  after  3000  s  for  250  m  resolution  and  lOth-order  polynomials,  (a) 
shows  the  total  domain  using  contour  values  between  —0.0015  and  0.003  with  a  contour  interval  of  0.0005  and  (b)  shows  the  profiles  along 
5000  m  height  for  all  five  models. 
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Table  2 

Case  1 :  Inertia-gravity  wave 


Model 

Wmax 

Wmin 

o' 

max 

Cin 

SE1 

2.698 

X 

1(T3 

-2.775 

X 

1(T3 

2.787 

X 

1(T3 

-1.518 

X 

1(T3 

SE2 

2.698 

X 

1(T3 

-2.774 

X 

1(T3 

2.787 

X 

1(T3 

-1.519 

X 

1(T3 

SE3 

2.698 

X 

1(T3 

-2.774 

X 

1(T3 

2.787 

X 

1(T3 

-1.519 

X 

1(T3 

DG2 

2.698 

X 

1(T3 

-2.774 

X 

1(T3 

2.787 

X 

1(T3 

-1.519 

X 

1(T3 

DG3 

2.698 

X 

1(T3 

-2.774 

X 

icr3 

2.787 

X 

1(T3 

-1.519 

X 

1(T3 

Comparison  of  the  five  models  studied  for  250  m  resolution  and  10th  order  polynomials  after  3000  s. 


there  is  a  mean  horizontal  flow  of  20  m/s  which  should  move  the  initial  perturbation  from  1 00,000  m  to 
160,000  m  in  3000  s;  this  is  in  fact  what  we  get  for  all  five  models. 

Skamarock  and  Klemp  [46]  give  an  analytic  solution  for  this  test  but,  unfortunately,  it  is  only  valid  for  the 
Boussinesq  equations,  which  while  useful  for  qualitative  comparisons,  cannot  be  used  to  compute  error  norms 
since  we  use  the  fully  compressible  equations.  We  use  the  same  contouring  interval  used  in  [46]  and  our  results 
look  quite  similar.  In  addition,  we  see  that  our  results  are  very  similar  to  those  of  the  recently  proposed  FV 
model  of  Ahmad  and  Lindeman  [1]  which  uses  equation  set  2.  They  use  a  smaller  contouring  interval  but  the 
potential  temperature  perturbations  of  our  models  have  similar  values.  Specifically,  their  values  are  in  the 
range  O'  e  [—1.49  x  10  3, 2.82  x  10  3]  whereas  ours  are  O'  e  [—1.51  x  10  3, 2.78  x  10  3]  which  we  show  in 
Table  2. 

Looking  at  Table  2  we  note  that  our  five  models  give  exactly  the  same  results  except  for  SE1.  We  do  not 
show  the  extrema  for  n'  and  u  because  they  are  identical  for  all  five  models;  for  n'  they  are 
(—1.238  x  10-6, 1.716  x  1CT6)  and  for  u  they  are  (—1.067  x  10  2, 1.069  x  10~2).  The  differences  in  wm in  and 
0'mm  for  SE1  are  certainly  quite  small  but  it  is  worth  noting  that  SE1  is  the  only  outlier.  The  results  for  this 
case  can  be  summarized  as  follows:  all  five  models  yield  very  good  results  and  sets  2  and  3  yield  identical  solu¬ 
tions,  with  set  1  being  the  only  one  that  yields  a  different  solution. 

6.2.  Case  2:  Rising  thermal  bubble 

Fig.  3  shows  the  potential  temperature  perturbation  contours  for  DG3  after  700  s  for  various  resolutions 
while  Fig.  4  shows  the  one-dimensional  profile  along  the  vertical  direction  for  x  =  500  m.  This  case  has  no 
analytic  solution  but  the  resulting  dynamics  are  sufficiently  simple  to  be  able  to  predict  its  proper  evolution. 
Fig.  4a  shows  that  at  a  20  m  resolution,  when  the  bubble  is  under-resolved  (see  Fig.  3a)  set  1  yields  a  different 
solution  from  sets  2  and  3.  Note  that  all  four  models  (SE2,  SE3,  DG2,  and  DG3)  give  the  exact  same  results 
and  SE1  being  the  only  outlier.  As  we  increase  the  resolution  to  10  m  (Fig.  4b)  and  to  5  m  (Fig.  4c),  the  DG 
models  (DG2  and  DG3)  yield  a  sharper  discontinuity  than  the  SE  models  (SE1,  SE2,  and  SE3).  The  sharper 
discontinuity  of  the  DG  models  can  also  be  seen  in  Table  3  where  the  solution  variables  are  tabulated  for  the 
5  m  resolution.  The  SE  models  yield  almost  identical  results  to  each  other  while  the  DG  models  yield  larger 
velocities  yet  smaller  potential  temperature  extrema  which  result  in  a  stronger  and  sharper  discontinuity. 
Finally,  at  a  3.5  m  resolution  (Fig.  4d)  all  the  models  converge  to  the  same  solution. 

The  results  for  this  test  can  be  summarized  as  follows:  set  1  should  not  be  chosen  for  under-resolved  flows; 
and  the  DG  models  are  preferable  to  the  SE  models. 

6.3.  Case  3:  Robert  smooth  bubble 

Fig.  5  shows  the  potential  temperature  perturbation  contours  for  DG3  after  700  and  800  s  for  a  5  m  reso¬ 
lution  while  Fig.  6  shows  the  one-dimensional  profile  along  the  vertical  direction  for  x  =  500  m  for  all  five 
models.  Fig.  5  shows  that  from  700  to  800  s,  the  dynamics  of  the  bubble  changes  from  very  smooth  to  turbu¬ 
lent.  Fig.  6a  shows  that  at  700  s,  the  SE  solutions  (SE1,  SE2,  and  SE3)  begin  to  separate  from  the  DG 
solutions  (DG2  and  DG3).  At  this  point,  it  is  not  known  which  solution  is  correct,  but  upon  looking  at 
Fig.  6b  we  see  that  at  800  s  the  SE  solutions  begin  to  oscillate  wildly  while  the  DG  solutions  maintain  a  sharp 
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Fig.  3.  Case  2:  Rising  Thermal  Bubble.  Potential  temperature  perturbation  after  700  s  for  the  following  resolutions:  (a)  20,  (b)  10,  (c)  5, 
and  (d)  3.5  m.  All  resolutions  use  10th  order  polynomials,  and  the  contour  values  are  from  —0.05  to  0.525  with  an  interval  of  0.025. 


discontinuity.  Table  4  confirms  that  the  DG  solutions  indeed  do  a  better  job  of  capturing  the  discontinuity 
since  it  can  be  observed  that  the  potential  temperature  (0')  extrema  are  smaller.  The  results  for  this  test  show 
that  for  strong  gradients,  the  DG  models  are  clearly  preferable  to  the  SE  models.  Note  that  we  have  not  used 
any  slope  limiter  or  additional  machinery  with  the  DG  method.  Both  the  SE  and  DG  models  share  the  exact 
same  machinery  such  as  the  spatial  filters.  This  test  highlights  the  advantages  that  the  DG  method  may  offer 
mesoscale  modeling.  In  future  work,  we  shall  include  slope  limiters  to  the  DG  method  in  order  to  improve  on 
the  handling  of  discontinuities  and  other  such  sharp  gradients. 

6.4.  Case  4:  Density  current 

Fig.  7  shows  the  potential  temperature  perturbation  contours  after  900  s  for  400,  200,  100,  and  50  m  res¬ 
olutions.  The  first  observation  is  that  even  at  the  very  coarse  resolution  of  400  m  (Fig.  7a),  two  of  the  three 
Kelvin-Helmholtz  rotors  are  clearly  visible.  At  a  200  m  resolution  (Fig.  7b)  the  second  rotor  is  much  better 
defined.  At  100  nr  (Fig.  7c),  all  three  rotors  are  clearly  visible.  Finally,  at  50  m  (Fig.  7d)  there  is  little  change 
in  the  structure  of  the  rotors,  the  only  discernible  difference  is  that  the  noise  above  the  rotors  has  disappeared. 

In  Fig.  8  we  plot  the  one-dimensional  profile  of  the  potential  temperature  perturbation  along  the  horizontal 
(x)  direction  at  a  height  of  z  =  1200  m  for  all  five  models  (on  the  left  panel)  and  the  profiles  for  five  different 
resolutions  (on  the  right  panel).  The  three  negative  wells  in  Fig.  8a  correspond  to  the  three  distinct  rotors 
clearly  visible  in  Fig.  7.  Fig.  8a  shows  that  the  results  for  the  five  models  split  into  two  distinct  groups.  In 
one  group,  the  three  models  SE1,  SE2,  and  DG2  all  give  the  same  results  and  in  the  second  group  the  models 
SE3  and  DG3  give  the  same  results.  The  differences  between  groups  one  and  two  are  quite  small  but  certainly 
not  insignificant.  This  difference,  however,  is  not  unexpected  since  SE3  and  DG3  use  equation  set  3  which  has 
the  true  Navier-Stokes  viscous  stress.  The  main  point  that  should  be  concluded  from  this  figure  is  that  the 
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Fig.  4.  Case  2:  Rising  Thermal  Bubble.  Potential  temperature  perturbation  after  700  s  for  the  following  resolutions:  (a)  20,  (b)  10,  (c)  5, 
and  (d)  3.5  m.  All  resolutions  use  10th  order  polynomials,  and  the  profiles  shown  are  along  the  vertical  (z)  at  x  =  500  m  for  all  five  models. 


Table  3 

Case  2:  Rising  thermal  bubble 


Model 

^max 

<nin 

^max 

^min 

Wmax 

Wmin 

O' 

max 

C, 

SE1 

9.364  x  1(T6 

-1.195  x  10-5 

2.073 

-2.073 

2.536 

-1.911 

0.570 

-0.098 

SE2 

9.364  x  10-6 

-1.196  x  ET5 

2.074 

-2.074 

2.536 

-1.911 

0.570 

-0.098 

SE3 

9.365  x  10-6 

-1.196  x  nr5 

2.074 

-2.074 

2.536 

-1.911 

0.570 

-0.098 

DG2 

9.355  x  10-6 

-1.195  x  1(T5 

2.081 

-2.081 

2.543 

-1.915 

0.538 

-0.093 

DG3 

9.355  x  1(T6 

-1.195  x  1(T5 

2.081 

-2.081 

2.543 

-1.915 

0.538 

-0.093 

Comparison  of  the  five  models  studied  for  5  m  resolution  and  10th  order  polynomials  after  700  s. 


equation  set  and  not  the  numerical  method  has  the  larger  impact  on  the  solution  for  this  test.  To  be  convinced 
of  this  one  only  needs  to  look  at  Table  5  which  we  will  discuss  in  a  moment.  Let  us  now  turn  to  Fig.  8b  to 
confirm  that  the  models  have  indeed  converged  at  100  m  resolution. 

Fig.  8b  shows  that  at  400  and  200  m  resolutions,  the  potential  temperature  profiles  are  still  changing.  How¬ 
ever,  looking  at  either  100,  50,  or  25  m  resolutions  one  can  see  that  all  these  curves  are  directly  on  top  of  each 
other;  which  confirms  that  our  models  have  converged  at  100  m  resolution.  Let  us  now  turn  to  the  results  of 
Table  5. 

Table  5  shows  that  the  results  for  SE2  and  DG2  are  in  fact  identical  and  the  results  for  SE3  and  DG3  are 
also  identical  to  each  other;  SE1  gives  different  results  to  SE2  and  SE3,  however,  these  differences  are  suffi¬ 
ciently  small  that  they  are  unnoticeable  in  the  one-dimensional  profile  plots. 

In  order  to  further  quantify  the  differences  between  the  five  models,  we  use  the  same  definition  for  the  front 
location  as  in  Straka  et  al.  [51],  that  is,  the  —1  °C  value  of  potential  temperature  perturbation.  In  Table  5  it  is 
evident  that  all  the  models  give  essentially  the  same  front  location;  for  set  1  the  front  occurs  earlier  than  for  the 
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Fig.  5.  Case  3:  Robert  Smooth  Bubble.  Potential  temperature  perturbation  after  (a)  700  and  (b)  800  s.  The  contour  values  shown  are  from 
—0.1  to  0.6  with  an  interval  of  0.05  using  10th  order  polynomials  and  a  5  m  resolution. 


Fig.  6.  Case  3:  Robert  Smooth  Bubble.  Potential  temperature  perturbation  profiles  after  (a)  700  and  (b)  800  s  for  all  five  models.  The 
profiles  shown  are  along  the  vertical  (r)  at  x  =  500  m  using  10th  order  polynomials  and  a  5  m  resolution. 


Table  4 

Case  3:  Robert  smooth  bubble 


Model 

n' 

max 

^max 

^min 

Wmax 

Wmin 

O’ 

max 

tflnin 

SE1 

5.994  x  1(T6 

-2.054  x  10-5 

2.471 

-2.471 

2.566 

-2.510 

0.625 

-0.233 

SE2 

5.993  x  10-6 

-2.065  x  10-5 

2.493 

-2.493 

2.566 

-2.529 

0.624 

-0.230 

SE3 

5.995  x  10-6 

-2.065  x  10-5 

2.493 

-2.493 

2.566 

-2.529 

0.624 

-0.230 

DG2 

5.979  x  10-6 

-2.025  x  10-5 

2.423 

-2.423 

2.553 

-2.335 

0.604 

-0.132 

DG3 

5.979  x  10-6 

-2.025  x  10~5 

2.423 

-2.423 

2.553 

-2.335 

0.604 

-0.132 

Comparison  of  the  five  models  studied  for  5  m  resolution  and  10th  order  polynomials  after  800  s. 


other  two  equation  sets.  Set  1  is  more  dissipative  as  indicated  by  the  smaller  potential  temperature  extrema. 
Equation  set  2  gives  less  dissipative  values  for  potential  temperature  and  the  front  location  thereby  occurs 
later.  Finally,  equation  set  3  yields  the  least  dissipative  solution  and  the  front  occurs  even  later  for  this  set. 
The  results  for  this  case  can  be  summarized  as  follows:  set  3  is  preferable  to  sets  1  and  2  -  the  choice  of 
SE  or  DG  is  not  as  important  as  the  equation  set. 

6.5.  Case  5:  Schcir  mountain 

Fig.  9  shows  that  the  numerical  (solid)  and  analytic  (dashed)  values  for  both  horizontal  and  vertical  veloc¬ 
ities  compare  quite  well.  In  Table  6,  we  compare  the  RMS  errors  for  all  four  variables  for  the  five  models. 
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a 


N 


c 


N 


4.800 

H  4.800 
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N  2.400 
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Fig.  7.  Case  4:  Density  Current.  Potential  temperature  perturbation  after  900  s  using  (a)  400,  (b)  200,  (c)  100,  and  (d)  50  m  resolution  with 
8th  order  polynomials.  The  contour  values  shown  are  between  —9  to  0  with  a  contour  interval  of  0.25. 


Fig.  8.  Case  4:  Density  Current.  Potential  temperature  perturbation  after  900  s  using  8th  order  polynomials  for  (a)  all  the  models  using 
25  m  resolution  and  (b)  DG3  at  five  different  resolutions  (400  m,  etc.).  Profiles  along  1200  m  height  are  illustrated. 


Table  5 

Case  4:  Density  current 


Model 

Front  location  (in  m) 

Tlf 

max 

"min 

O' 

max 

SE1 

14,736 

8.431  x  1(T4 

-1.662  x  1(T2 

0.00059 

-8.83 

SE2 

14,767 

8.428  x  1(T4 

-1.657  x  1(T2 

0.00012 

-8.90 

SE3 

14,789 

8.361  x  1(T4 

-1.690  x  10~2 

0.00331 

-9.08 

DG2 

14,767 

8.428  x  1(T4 

-1.657  x  10~2 

0.00012 

-8.90 

DG3 

14,789 

8.361  x  1(T4 

-1.690  x  10“2 

0.00331 

-9.08 

Comparison  of  the  five  models  studied  for  25  m  resolution  and  8th  order  polynomials  after  900  s. 


Clearly,  all  five  models  give  essentially  identical  results;  the  differences  between  the  models  are  quite  small. 
Nonetheless,  the  results  for  SE1  are  slightly  better  than  those  for  SE2  and  SE3.  Looking  at  the  rest  of  the  mod¬ 
els  we  see  that  DG2  and  DG3  give  better  results  than  any  of  the  SE  models.  For  this  test,  the  conclusion  is  that 
one  should  use  DG  over  SE. 

6.6.  Case  6:  Linear  hydrostatic  mountain 

Fig.  10  shows  that  the  numerical  (solid)  and  analytic  (dashed)  values  for  both  horizontal  and  vertical  veloc¬ 
ities  compare  quite  well.  Clearly,  the  vertical  velocity  is  captured  better  than  the  horizontal  velocity  as  can  be 
seen  in  Table  7  where  we  see  that  the  RMS  errors  differ  by  almost  two  orders  of  magnitude.  This  is  not 
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Fig.  9.  Case  5:  Schar  Mountain.  The  numerical  solution  (solid  line)  and  analytic  solution  (dashed  line)  after  10  h  for  250  m  (in  x)  and 
210  m  (in  r)  resolution  and  10th  order  polynomials.  For  the  horizontal  velocity  (a)  the  contour  values  are  between  —2  to  +2  with  a  contour 
interval  of  0.2.  For  the  vertical  velocity  (b)  the  contour  values  are  between  —2  to  +2  with  a  contour  interval  of  0.05. 


Table  6 

Case  5:  Schar  mountain 


Variable 

SEl 

SE2 

SE3 

DG2 

DG3 

71 

8.20  x  10~6 

8.29  x  10~6 

8.27  x  10~6 

7.37  x  10~6 

7.36  x  10~6 

U 

2.26  x  10“' 

2.26  x  10~‘ 

2.26  x  10  1 

1.94  x  10*1 

1.94  x  10  1 

W 

7.58  x  10~2 

7.67  x  10~2 

7.66  x  10~2 

7.51  x  10-2 

7.51  x  10~2 

e 

6.88  x  10~2 

6.80  x  10~2 

6.78  x  10~2 

5.85  x  10~2 

5.84  x  10~2 

Root-mean-square  errors  for  all  four  variables  after  10  h  for  250  m  (in  x)  and  210  m  (in  r)  resolution  and  10th  order  polynomials  for  all 
five  models. 


Fig.  10.  Case  6:  Linear  Flydrostatic  Mountain.  The  numerical  solution  (solid  line)  and  analytic  solution  (dashed  line)  after  10  h  for  1200  m 
(in  x)  and  240  m  (in  r)  resolution  and  10th  order  polynomials.  For  the  horizontal  velocity  (a)  the  contour  values  are  between  —0.025  to 
+0.025  with  a  contour  interval  of  0.005.  For  the  vertical  velocity  (b)  the  contour  values  are  between  —0.005  to  +0.005  with  a  contour 
interval  of  0.0005. 


surprising  for  this  case  since  it  is  in  the  hydrostatic  regime  and,  thereby,  the  vertical  acceleration  does  not  play 
a  dominant  role. 

There  are  a  number  of  interesting  observations  that  can  be  discerned  from  Table  7.  The  first  is  that  all 
five  models  give  relatively  similar  results.  The  second  observation  is  that  the  DG  models  give  slightly  bet¬ 
ter  results  than  the  SE  models;  one  can  clearly  see  this  by  comparing  the  RMS  errors  for  n1,  u,  and  w 
(except  for  6').  The  third  observation  is  that  equation  sets  2  and  3  give  better  results  than  equation  set 
1.  Note  that  SE1,  SE2,  and  SE3  all  use  spectral  elements  -  the  only  difference  between  these  three  models 
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Table  7 

Case  6:  Linear  hydrostatic  mountain 


Variable 

SEl 

SE2 

SE3 

DG2 

DG3 

71 

1.21  x  1(T7 

1.21  x  1(T7 

1.21  x  1(T7 

1.20  x  10~7 

1.20  x  10~7 

U 

2.28  x  1(T3 

2.24  x  1(T3 

2.24  x  1(T3 

2.23  x  10-3 

2.23  x  10~3 

W 

5.98  x  10-5 

4.31  x  1(T5 

4.31  x  1(T5 

3.20  x  10~5 

3.20  x  10~5 

0 

1.67  x  1CT3 

1.64  x  1(T3 

1.64  x  1CT3 

1.65  x  10“3 

1.65  x  10~3 

Root-mean-square  errors  for  all  four  variables  for  1200  m  (in  .v)  and  240  m  (in  z)  resolution  and  10th  order  polynomials  after  10  h  for  all 
five  models. 


is  the  equation  sets  used.  Finally,  the  last  observation  is  that  equation  sets  2  and  3  yield  identical  results 
when  using  the  same  numerical  method. 

In  Fig.  1 1  we  plot  the  normalized  momentum  flux  for  DG3  at  various  times  in  the  integration  (Fig.  1  la)  and 
for  all  five  models  after  10  h  (Fig.  lib).  Fig.  11a  shows  that  the  simulations  have  reached  steady-state  after 
10  h.  Fig.  lib  shows  that  the  normalized  momentum  flux  values  are  almost  identical  for  all  five  models 
and  are  quite  good,  that  is,  the  values  are  between  0.95  and  1.01  -  the  best  one  could  do  is  to  have  values 
of  unity  which  we  approach  quite  quickly.  The  results  of  this  test  case  can  be  summarized  as  follows:  all  things 
being  equal,  one  should  use  sets  2  or  3  instead  of  set  1  and  DG  is  slightly  better  than  SE. 

6. 7.  Case  7:  Linear  nonhydrostatic  mountain 

Fig.  12  shows  the  horizontal  and  vertical  velocity  contours  for  DG3  after  5  h.  Note  how  different  the  ver¬ 
tical  velocity  fields  look  between  the  nonhydrostatic  (Fig.  12)  and  hydrostatic  (Fig.  10)  cases.  In  Fig.  12  we 
compare  the  numerical  solution  (solid  lines)  with  the  analytic  solution  (dashed  lines).  The  simulated  horizontal 
velocity  (Fig.  12a)  compares  quite  well  with  the  analytic  solution  near  the  mountain  but  loses  accuracy  away 
from  the  mountain.  In  contrast,  the  simulated  vertical  velocity  (Fig.  12b)  compares  extremely  well  with  the 
analytic  solution  everywhere. 

In  Fig.  13  we  plot  the  normalized  momentum  flux  for  DG3  at  various  times  in  the  integration  (Fig.  13a)  and 
for  all  five  models  after  5  h  (Fig.  13b).  Fig.  13a  shows  that  the  simulations  have  reached  steady-state  after  5  h; 
this  is  clearly  seen  by  the  convergence  of  the  momentum  fluxes  for  times  3,  4,  and  5  h.  Fig.  13b  shows  that  all 
five  models  give  very  similar  momentum  flux  values.  In  fact,  the  values  are  very  good  since  they  remain 
between  0.95  and  1. 

In  Table  8  we  list  the  RMS  errors  for  all  four  variables  of  the  five  models.  These  results  show  that  all  five 
models,  once  again,  yield  very  similar  results.  Although  the  differences  in  the  RMS  errors  are  quite  small  it  can 
be  seen  that:  SE2  and  SE3  give  better  results  than  SE1;  and  DG2  and  DG3  give  better  results  than  the  SE 


Normalized  Momentum  Flux 


Normalized  Momentum  Flux 


Fig.  11.  Case  6:  Linear  Hydrostatic  Mountain.  Normalized  momentum  flux  for  1200  m  (in  x)  and  240  m  (in  z)  resolution  and  10th  order 
polynomials  for  (a)  DG3  at  times  2,  4,  6,  8,  and  10  h,  and  (b)  for  all  five  models  at  10  h. 
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a  12.000 
10.000 

8.000 

n  6.000 

4.000 

2.000 

0 

60.000  75.000  90.000  105.000  60.000  75.000  90.000  105.000 

x  x 

Fig.  12.  Case  7:  Linear  Nonhydrostatic  Mountain.  The  numerical  solution  (solid  line)  and  analytic  solution  (dashed  line)  after  5  h  for 
360  m  (in  a)  and  300  m  (in  z)  resolution  and  10th  order  polynomials.  For  the  horizontal  velocity  (a)  the  contour  values  are  between  —0.025 
to  +0.025  with  a  contour  interval  of  0.0025.  For  the  vertical  velocity  (b)  the  contour  values  are  between  —0.005  to  +0.005  with  a  contour 
interval  of  0.0005. 


Fig.  13.  Case  7:  Linear  Nonhydrostatic  Mountain.  Normalized  momentum  flux  for  360  m  (in  x )  and  300  m  (in  z)  resolution  and  10th 
order  polynomials  for  (a)  DG3  at  times  1,  2,  3,  4,  and  5  h,  and  (b)  for  all  five  models  at  5  h. 


Table  8 

Case  7:  Linear  nonhydrostatic  mountain 


Variable 

SEl 

SE2 

SE3 

DG2 

DG3 

7C 

1.69  x  10~8 

1.62  x  10~8 

1.62  x  10~8 

1.65  x  10~8 

1.61  x  10~8 

U 

4.85  x  10~4 

4.84  x  10“4 

4.84  x  10~4 

4.79  x  10~4 

4.80  x  10~4 

w 

9.39  x  10-5 

9.39  x  10'5 

9.39  x  10-5 

6.57  x  10-5 

8.94  x  10~5 

e 

1.69  x  10~4 

1.69  x  10“4 

1.69  x  10~4 

1.67  x  10~4 

1.68  x  10~4 

Root-mean-square  errors  for  all  four  variables  for  360  m  (in  x)  and  300  m  (in  z)  resolution  and  10th  order  polynomials  after  10  h  for  all 
five  models. 


models.  The  results  for  this  test  can  be  summarized  as  follows:  all  models  performed  quite  similarly,  however, 
sets  2  and  3  are  preferred  over  set  1  and  DG  is  preferred  over  SE. 

6.8.  Conservation  and  efficiency 

A  few  words  regarding  the  conservation  properties  and  efficiency  of  the  models  are  in  order.  We  have  dis¬ 
cussed  the  three  equation  sets  and  have  mentioned  that  only  sets  2  and  3  are  formally  conservative.  However, 
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Fig.  14.  Conservation  of  (a)  mass  and  (b)  energy  for  all  five  models  using  the  strong  form.  The  results  shown  are  for  the  rising  thermal 
bubble  using  a  20  m  resolution  with  10th  order  polynomials. 


the  question  that  needs  to  be  addressed  is  what  exactly  are  these  equations  conserving.  For  equation  set  3  the 
answer  is  obvious  since  the  governing  equations  consist  of  three  conservation  laws  for  mass,  momentum,  and 
total  energy.  For  set  2,  the  conserved  quantities  are  mass,  momentum,  and  density  potential  temperature.  For 
set  1,  none  of  the  variables  should  be  conserved.  Fig.  14  shows  the  loss  of  mass  (left  panel)  and  energy  (right 
panel)  conservation  for  all  five  models  using  the  strong  form  of  the  equations  (Eqs.  (23)  and  (32)).  As 
expected,  Fig.  14a  shows  that  the  only  model  that  does  not  conserve  mass  is  set  1  (SE1);  it  should  be  under¬ 
stood  that  this  is  not  a  problem  with  the  numerical  method  but  rather  with  the  governing  equation  -  this  form 
cannot  formally  conserve  any  quantity.  Fig.  14b  shows  that  the  only  models  which  conserve  energy  are  SE3 
and  DG3,  that  is,  those  that  use  set  3.  On  the  other  hand,  SEE  SE2,  and  DG2  do  not  conserve  energy.  We 
know  already  that  set  1  will  not  conserve  any  quantity,  however,  set  2  does  not  conserve  total  energy  but  it 
does  in  fact  conserve  density  potential  temperature  (not  shown).  Potential  temperature  is  directly  related  to 
entropy  and  therefore  it  makes  more  sense  to  conserve  energy  than  to  impose  the  model  to  be  isentropic.  Per¬ 
haps  once  physical  parameterization  (in  the  form  of  source  functions  in  momentum  and  energy)  is  added  to 
the  equations  to  represent  sub-grid  scale  processes,  the  debate  whether  conserving  total  energy  or  entropy 
becomes  moot.  We  plan  to  revisit  this  issue  in  future  work  where  we  shall  experiment  with  simple  physical 
parameterizations. 

Note  that  while  the  energy  conservation  given  in  Fig.  14b  is  quite  good,  it  is  not  conservative  up  to  machine 
precision.  This  is  due  to  the  inexact  integration  we  use  for  both  the  SE  and  DG  methods.  Flowever,  if  the  weak 
form  of  the  equations  are  used  (Eqs.  (25)  and  (31)),  then  SE3  and  DG3  are  conservative  up  to  machine  pre¬ 
cision  (10~16)  which  we  show  in  Fig.  15b. 


Fig.  15.  Conservation  of  (a)  mass  and  (b)  energy  for  all  five  models  using  the  weak  form.  The  results  shown  are  for  the  rising  thermal 
bubble  using  a  20  m  resolution  with  10th  order  polynomials. 
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Fig.  16.  Efficiency  of  the  (a)  strong  and  (b)  weak  forms  of  the  equations  for  all  five  models.  The  results  shown  are  for  the  rising  thermal 
bubble  using  various  number  of  elements  in  the  horizontal  (with  equal  number  in  the  vertical)  using  10th  order  polynomials  with  an 
integration  time  of  700  s. 


While  we  have  not  devoted  a  lengthy  discussion  to  the  efficiency  comparisons  of  the  five  models,  we  should 
mention  that  the  DG  models  are,  in  general,  more  computationally  expensive  than  the  SE  models.  The  reason 
for  the  heavier  costs  incurred  by  DG  methods  is  due  to  the  fact  that  the  DG  method  requires  the  computation 
of  area  integrals  typical  of  SE  methods  as  well  as  boundary  integrals  typical  of  FV  methods.  The  cost  of  area 
integrals  is  0(NeN 3)  where  Ne  are  the  number  of  elements  and  N  is  the  order  of  the  polynomial.  The  cost  of 
boundary  integrals  is  0(NsN)  where  Ns  is  the  number  of  sides/edges.  For  viscous  problems,  additional  floating 
point  operations  of  O (NeN*)  are  required  by  the  local  discontinuous  Galerkin  (LDG)  method.  To  give  the 
reader  a  sense  of  these  costs,  in  Fig.  16  we  plot  the  wallclock  time  required  by  all  five  models  to  integrate 
700  s  for  the  rising  thermal  bubble  for  various  resolutions.  Recall  that  the  rising  thermal  bubble  (case  2)  is 
solved  on  a  square  domain  with  equal  resolution  in  the  horizontal  (x)  and  vertical  (z)  directions.  In 
Fig.  16a  we  plot  the  wallclock  time  for  all  five  models  using  the  strong  form  of  the  equations  which  are 
not  guaranteed  to  conserve  exactly  even  for  sets  2  and  3.  We  should  point  out,  however,  that  this  is  the  worst 
case  scenario  for  the  DG  method  -  meaning  that  for  these  timings,  we  have  used  the  maximum  time-steps 
allowed  by  the  SE  models  which  are  50%  larger  than  those  allowed  by  the  DG  models. 

If  we  impose  conservation,  then  the  equations  have  to  be  solved  in  weak  form  and  the  maximum  allowable 
time-step  for  the  SE  method  is  less  than  or  equal  to  that  for  the  DG  method;  this  slightly  larger  time-step 
allows  the  DG  method  (although  it  requires  more  floating  point  operations)  to  compete  with  the  SE  method. 
Comparing  the  weak  form  solutions,  Fig.  16b  shows  that  the  DG  models  are  now  as  efficient  as  the  SE  models; 
note  that  we  have  also  plotted  the  SE1  solution  although  this  form  is  not  conservative  even  in  the  weak  form. 
Another  observation  regarding  the  efficiency  of  the  three  equation  sets  is  that  set  3  is  more  efficient  to  compute 
than  set  2,  and  set  1  is  more  efficient  than  sets  2  and  3.  The  reason  has  to  do  with  the  cost  of  computing  the 
pressure.  For  set  1,  the  pressure  (Exner  in  this  case)  is  in  fact  a  prognostic  variable  whereas  for  sets  2  and  3  it  is 
diagnostic  and  hence  has  to  be  computed  from  the  prognostic  variables.  For  set  2,  the  pressure  is  a  function  of 
the  prognostic  variables  raised  to  an  exponent  (see  Eq.  (5))  whereas  for  set  3,  the  pressure  is  a  linear  function 
of  the  prognostic  variables  (see  Eq.  (8))  and  therefore  cheaper  to  compute.  This  brief  discussion  on  the  cost  of 
the  SE  and  DG  methods  is  completely  independent  from  the  time-integrators  used.  In  this  work  we  only  dis¬ 
cuss  explicit  methods  but  in  a  separate  paper  we  present  semi-implicit  results  [41], 

7.  Conclusions 

In  this  paper  we  studied  three  forms  of  the  Euler  equations  for  possible  use  in  next-generation  nonhydro¬ 
static  models.  Set  1  is  not  in  conservation  form;  set  2  can  be  written  in  conservation  form  only  for  inviscid 
flow;  while  set  3  can  be  written  in  conservation  form  for  either  inviscid  or  viscous  flow  applications  and  is 
in  fact  the  compressible  Navier-Stokes  (NS)  equations.  We  compared  the  results  of  the  three  equation  sets 


F.X.  Giraldo,  M.  Restelli  I  Journal  of  Computational  Physics  227  (2008)  3849-3877 


3875 


using  two  numerical  methods:  the  spectral  element  (SE)  and  discontinuous  Galerkin  (DG)  methods;  both  are 
local  high-order  methods  but  the  DG  method  is  a  fully  conservative  method  (both  global  and  local)  while  the 
SE  method  is  only  globally  conservative.  We  tested  the  five  models  using  seven  test  cases:  three  of  which 
involved  mountain  waves.  The  live  models  gave  similar  results  and  compare  well  with  models  found  in  the 
literature.  Equation  set  1,  although  not  in  conservation  form,  performed  better  than  expected  although  it 
tended  to  yield  slightly  different  results  than  sets  2  and  3.  The  differences  between  the  five  models  were  quite 
small  but  the  main  result  found  is  that  the  differences  in  the  simulation  were  most  often  attributed  to  the  equa¬ 
tion  set  used;  the  numerical  method  also  played  a  role  in  the  performance  of  the  models.  Set  3  (fully  compress¬ 
ible  NS  equations)  yielded  solutions  that  were  less  dissipative  (for  the  density  current).  The  comparison  of  the 
SE  and  DG  methods  showed  that  both  methods  yield  similar  solutions  for  the  types  of  problems  studied;  how¬ 
ever,  for  the  bubble  and  mountain  tests,  the  DG  method  performed  better.  It  is  expected  that  for  more  strin¬ 
gent  tests  (i.e.  those  having  very  steep  gradients  such  as  those  produced  by  the  sub-grid  scale  parameterization) 
the  DG  method  will  prevail.  In  future  work  we  will  explore  high-order  non-reflecting  boundary  conditions  and 
add  sub-grid  scale  parameterization  to  study  the  sensitivity  of  these  methods  to  the  sharp  local  gradients  it 
produces.  Although  it  is  too  early  to  conclude  which  equation  set  and  method  will  work  best  with  the  param¬ 
eterization,  we  expect  the  DG  method  to  handle  the  sub-grid  scale  processes  better  due  to  its  ability  to  better 
handle  discontinuities  (as  in  the  two  bubble  problems).  The  exact  conservation  properties  of  set  3  in  its  weak 
form  motivates  us  to  pursue  this  equation  set  in  our  future  studies. 
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